home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / fp_compute.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-27  |  5.8 KB  |  214 lines
  1. /*
  2. ### Driver for computing periodic orbits of a given period ###
  3.  
  4. Note: All computation in primary coordinates
  5. */
  6.  
  7. #include "../modellib/class_kaos_def.h"
  8.  
  9. fp_compute()
  10. {
  11.     int i,j,k,kc,n,color,fp_type,s_type,ndid,duplicate_found;
  12.     static int ii=1;
  13.     double fabs(),sum,xerr,time,*x1,*x2,*xt,*xo,**evec,**eval,*dvector(),**dmatrix();
  14.     extern int stop,var_dim,ir,symbol_size,fp_show_toggle,fp_display_option,fp_go_option,fp_algorithm,polar_coord,sqzex_maxsq;
  15.     extern int cur_color,polar_coord,mapping_on,model,enable_period;
  16.         extern int n_stored_fp,mdim_fp,mdim_fp_init,mdim_fp_crit,mdim_fp_step,n_mc,sqzex_maxsq,fi_maxsq;
  17.     extern double mc_eps,sqzex_eps,sqzex_epsf,*period_len,*win_var_i,*all_max,*all_min,*param;
  18.         extern char string[],string2[];
  19.     extern int *fp_period;
  20.     extern double **fp_x,***fp_eval,***fp_evec,*fp_xerr;
  21.     extern int (*f_p)();
  22.  
  23.         
  24.     stop = 0;
  25.     n=var_dim;
  26.     x1 = dvector(0,n-1);
  27.     x2 = dvector(0,n-1);
  28.     xt = dvector(0,n-1);
  29.     xo = dvector(0,n-1);
  30.     eval = dmatrix(0,n-1,0,1);
  31.     evec = dmatrix(0,n-1,0,n-1);
  32.  
  33.     /* ir: period of an orbit */
  34.         if(fp_go_option==0){
  35.                 n_mc=1;        /* # of Newton try = 1 using a Mouse input*/
  36.         }
  37.  
  38.     if(!mapping_on) {
  39.         if(ir !=0) {
  40.             system_mess_proc(1,"Period>0 not implemented yet! Set period=0.");
  41.             ir = 0;
  42.         }
  43.     }
  44.     else {
  45.         if(ir ==0){
  46.             system_mess_proc(1,"Period=0 is not acceptable. Change period in the periodic orbit window (Period>0).");
  47.             goto done;
  48.  
  49.         }
  50.     }
  51.     sprintf(string,"Period=%d",ir);
  52.     system_mess_proc(0,string);
  53.  
  54.         /* Monte Carlo computation of all fixed points */
  55.  
  56.     srandom(ii++);
  57.     if(fp_algorithm==0)
  58.         system_mess_proc(0,"Starting Monte-Carlo Newton search...");
  59.     else
  60.         system_mess_proc(0,"Starting Monte-Carlo Secant search...");
  61.     for(kc=0;kc<n_mc;kc++){
  62.         /* get initial conditions from windows */
  63.         if(fp_go_option == 0){
  64.             for(i=0;i<n;i++) x1[i] = win_var_i[i];
  65.         }
  66.         else {
  67.             for(i=0;i<n;i++){
  68.                 x1[i] = (double) (random() & 0777777) / 0777777;
  69.                 x1[i] = all_min[i] + x1[i] * (all_max[i] - all_min[i]);
  70.             }
  71.         }
  72.         /* all computation in enclidean coords */
  73.         from_window_variables(xo,x1,polar_coord);
  74.  
  75.         if(fp_algorithm==0){
  76.             mnewt(xo,sqzex_maxsq,n,ir,sqzex_eps,sqzex_epsf,&xerr,&ndid);
  77.         }
  78.         else {
  79.             msecant(xo,sqzex_maxsq,n,ir,sqzex_eps,sqzex_epsf,&xerr,&ndid);
  80.         }
  81.         if(stop==1){
  82.             goto done;
  83.         }
  84.         if(xerr < sqzex_eps) {
  85.             for(i=0;i<n;i++) x1[i] = xo[i];
  86.             /* eliminate duplicates */
  87.             duplicate_found=0;
  88.             for(k=0;k<n_stored_fp;k++) {
  89.                 for(i=0;i<n;i++)xt[i] = x1[i] - fp_x[i][k];
  90.                 if(enable_period)
  91.                     dist_periodic(x2,xt,n);
  92.                 else
  93.                     for(i=0;i<n;i++)x2[i] = xt[i];
  94.                     
  95.                 sum =0;
  96.                 for(i=0;i<n;i++) sum += fabs(x2[i]);
  97.                 if(sum < mc_eps) {
  98.                     duplicate_found=1;
  99.                     break;
  100.                 }
  101.                 /* test for other members of a periodic orbit*/
  102.                 if(mapping_on) {
  103.                     for(i=1;i<ir;i++){
  104.                         (int) f_p(x2,1,x1,param,time,n);
  105.                         for(j=0;j<n;j++) x1[j] = x2[j];
  106.                         for(j=0;j<n;j++)xt[j] = x1[j] - fp_x[j][k];
  107.                         if(enable_period)
  108.                             dist_periodic(x2,xt,n);
  109.                         else
  110.                             for(i=0;i<n;i++)x2[i] = xt[i];
  111.                         sum =0;
  112.                         for(j=0;j<n;j++) sum += fabs(x2[j]);
  113.                         if(sum < mc_eps) {
  114.                             duplicate_found=1;
  115.                             break;
  116.                         }
  117.                     }
  118.                     if(duplicate_found)
  119.                         break;
  120.                 }
  121.                 else {
  122.                     /* ode version need to be installed */
  123.                 }
  124.             }
  125.  
  126.             if(duplicate_found){
  127.                     /* Print distinguished fixed points */
  128.                 printf("Dup PO: Per=%d NDid=%d X=(",ir,ndid); 
  129.                 for(i=0;i<n;i++) printf("%g,",xo[i]);
  130.                 printf("),Err=%g\n",xerr); 
  131.             }
  132.             else {
  133.                     /* Print distinguished fixed points */
  134.                 printf("New PO: Per=%d NDid=%d X=(",ir,ndid); 
  135.                 for(i=0;i<n;i++) printf("%g,",xo[i]);
  136.                 printf("),Err=%g\n",xerr); 
  137.                 /* record data */
  138.                 for(i=0;i<n;i++) fp_x[i][n_stored_fp]=xo[i];
  139.                 fp_xerr[n_stored_fp] = xerr;
  140.                 fp_period[n_stored_fp] = ir;
  141.                 if(n_stored_fp >= mdim_fp-1) {
  142.                     if(mdim_fp < mdim_fp_crit)
  143.                         mdim_fp *= 2;
  144.                     else
  145.                         mdim_fp += mdim_fp_step;
  146.                     realloc_fp_data();
  147.                 }
  148.  
  149.                 (void) cycle_color(&cur_color);
  150.                 (void) draw_record_pwf(xo,cur_color,4,5,1,0);
  151.                 (void) draw_record_pwf(xo,cur_color,4,5,1,0);
  152.                 if(fp_display_option==1){
  153.                     /* record_on = 0: do not record
  154.                     copies of a periodic orbit */
  155.                     if(ir>1) {
  156.                         for(i=0;i<n;i++) x1[i] = xo[i];
  157.                         (void) draw_record_other_pwf(x1,ir,cur_color,4,5,1,0);
  158.                         for(i=0;i<n;i++) x1[i] = xo[i];
  159.                         (void) draw_record_other_pwf(x1,ir,cur_color,4,5,1,0);
  160.                     }
  161.                 }
  162.  
  163.                 /* get eigen info */
  164.                 fp_get_evinfo(&fp_type,eval,evec,xo,ir,n);
  165.                 /* record all eigenvalues */
  166.                 for(i=0;i<n;i++){
  167.                     for(j=0;j<2;j++) fp_eval[i][j][n_stored_fp] = eval[i][j];
  168.                 }
  169.  
  170.                 /* record only saddle eigenvectors only (temporary)*/
  171.                 if(fp_type==0){
  172.                     for(i=0;i<n;i++){
  173.                         for(j=0;j<n;j++) fp_evec[i][j][n_stored_fp] = evec[i][j];
  174.                     }
  175.                 }
  176.                 /* get label, color and symbol_type for each
  177.                 fixed point type */
  178.                 fp_get_attributes(string,&color,&s_type,fp_type);
  179.                 (void) draw_record_pwf(xo,color,s_type,symbol_size,1,0);
  180.                 if(fp_display_option==1){
  181.                     /* record_on = 0: do not record
  182.                     copies of a periodic orbit */
  183.                     if(ir>1) {
  184.                         for(i=0;i<n;i++) x1[i] = xo[i];
  185.                         (void) draw_record_other_pwf(x1,ir,color,s_type,symbol_size,1,0);
  186.                     }
  187.                 }
  188.  
  189.                 printf("   Type=%s,EVal=",string);
  190.                  for(j=0;j<n;j++) printf(" (%g, %g)",fp_eval[j][0][n_stored_fp],fp_eval[j][1][n_stored_fp]);
  191.                 printf("\n");
  192.                 /* reset panel_item */
  193.                 for(j=0;j<n;j++) x1[j]=fp_x[j][n_stored_fp];
  194.                  to_window_variables(win_var_i,x1,polar_coord);
  195.                 n_stored_fp++;
  196.             }
  197.  
  198.         }
  199.         else {
  200.             printf("Bad PO: Per=%d NDid=%d,Err=%g\n",ir,ndid,xerr);
  201.         }
  202.         cur_color++;
  203.     }
  204.  
  205. done:
  206.     system_mess_proc(0,"Done!");
  207.     (void) free_dvector(x1,0,n-1);
  208.     (void) free_dvector(x2,0,n-1);
  209.     (void) free_dvector(xt,0,n-1);
  210.     (void) free_dvector(xo,0,n-1);
  211.     (void) free_dmatrix(eval,0,n-1,0,1);
  212.     (void) free_dmatrix(evec,0,n-1,0,n-1);
  213. }
  214.